GLOE-Seq analysis

library(wigglescout)
library(ggplot2)

RStudio Community is a great place to get help: https://community.rstudio.com/c/tidyverse

library(reshape2)
library(dplyr)

Attaching package: ‘dplyr’

The following objects are masked from ‘package:stats’:

    filter, lag

The following objects are masked from ‘package:base’:

    intersect, setdiff, setequal, union

library(ggpubr)

Registered S3 method overwritten by 'data.table':
  method           from
  print.data.table     

library(GenomicRanges)

Loading required package: stats4
Loading required package: BiocGenerics
Loading required package: parallel

Attaching package: ‘BiocGenerics’

The following objects are masked from ‘package:parallel’:

    clusterApply, clusterApplyLB, clusterCall, clusterEvalQ, clusterExport, clusterMap, parApply,
    parCapply, parLapply, parLapplyLB, parRapply, parSapply, parSapplyLB

The following objects are masked from ‘package:dplyr’:

    combine, intersect, setdiff, union

The following objects are masked from ‘package:stats’:

    IQR, mad, sd, var, xtabs

The following objects are masked from ‘package:base’:

    anyDuplicated, append, as.data.frame, basename, cbind, colnames, dirname, do.call, duplicated, eval,
    evalq, Filter, Find, get, grep, grepl, intersect, is.unsorted, lapply, Map, mapply, match, mget, order,
    paste, pmax, pmax.int, pmin, pmin.int, Position, rank, rbind, Reduce, rownames, sapply, setdiff, sort,
    table, tapply, union, unique, unsplit, which.max, which.min

Loading required package: S4Vectors

Attaching package: ‘S4Vectors’

The following objects are masked from ‘package:dplyr’:

    first, rename

The following object is masked from ‘package:base’:

    expand.grid

Loading required package: IRanges

Attaching package: ‘IRanges’

The following objects are masked from ‘package:dplyr’:

    collapse, desc, slice

Loading required package: GenomeInfoDb

source('../lib/lib_embed_data.v1.0.R')

Loading required package: xfun

Attaching package: ‘xfun’

The following objects are masked from ‘package:base’:

    attr, isFALSE

data_path <- '/Volumes/DATA/DATA/GLOE/GLOE/replicates/'

replicates <- c("WT-UT-1-1","WT-UT-1-2","WT-UT-1-3","WT-20-2-1","WT-20-2-2","WT-20-2-3","WT-24hr-3-1","WT-24hr-3-2","WT-24hr-3-3","KO-UT-4-1","KO-UT-4-2","KO-UT-4-3","KO-20-5-1","KO-20-5-2","KO-20-5-3","KO-24hr-6-1","KO-24hr-6-2","KO-24hr-6-3")
samples <- c("WT-UT-1","WT-20-2","WT-24hr-3","KO-UT-4","KO-20-5","KO-24hr-6")

build_filenames <- function(path, cond) {
  files <- paste0(path,"/GL-",cond,".hg19.bs50.bw")
  labels <- paste0("GLOE ", cond)
  return(data.frame(files=files,labels=labels))
}

build_filenames_plus <- function(path, cond) {
  files <- paste0(path,"/GL-",cond,".hg19.R1plus.bw")
  labels <- paste0("GLOE (+) ", cond)
  return(data.frame(files=files,labels=labels))
}

build_filenames_minus <- function(path, cond) {
  files <- paste0(path,"//GL-",cond,".hg19.R1minus.bw")
  labels <- paste0("GLOE (-) ", cond)
  return(data.frame(files=files,labels=labels))
}

Published peaksets (Cistrome)

Here, we are comparing profiles accross known accessible, DNAase hypersensitive regions, and CTCF peaks Note: unlike in previous plots, we clearly see now a correlation of CAD activity and accessible chromatin. Plotting average profiles on scarce read data can be thrown of by artifact regions. These artifacts are avoided here by 1) using a blacklist to remove known repetitive/artifact-prone regions and 2) removing outliers (top 0.1%) of the signal before averaging.

fn <- build_filenames(data_path,samples)

p1 <- plot_bw_profile(bwfiles = fn$files,labels = fn$labels,mode = "center",upstream = 1000, downstream = 1000, loci = "../genome/HCT116_ATACseq.hg19.bed", bin_size = 10 ,remove_top = 0.001, verbose=F)

p2 <- plot_bw_profile(bwfiles = fn$files,labels = fn$labels,mode = "center",upstream = 1000, downstream = 1000, loci = "../genome/HCT116_CTCF.hg19.bed", bin_size = 10 ,remove_top = 0.001, verbose=F)

p3 <- plot_bw_profile(bwfiles = fn$files,labels = fn$labels,mode = "center",upstream = 1000, downstream = 1000, loci = "../genome/HCT116_DHS.hg19.bed", bin_size = 10 ,remove_top = 0.001, verbose=F)

ggarrange(p1, p2, p3, ncol = 3, nrow = 1)

Pileup regions = ‘hypersensitive breakpoints’

Given the overall low coverage/incidence of SSB in the genome, a coincidence of several SSBs in the same location (or even exact same base pair) is highly unlikely to occur by chance. Noticing that around CTCF sites, uniqe GLOE-Seq reads piled up neatly (some with exact same starting base), I developed a pipeline that looks for such pileups. There are ~40 of these regions identified, which I would call ‘hypersensitive’. There may be many more but it is not possible to call them confidently.

fn <- build_filenames(data_path,samples)

p1 <- plot_bw_profile(bwfiles = fn$files,labels = fn$labels,mode = "center",upstream = 1000, downstream = 1000, loci = "../data/pileup/GL-WT-24hr-3.pile_up3_s10_u.uniq.bed", bin_size = 10, verbose = F) + scale_y_continuous(limits=c(0,45))

p3 <- plot_bw_profile(bwfiles = fn$files,labels = fn$labels,mode = "center",upstream = 1000, downstream = 1000, loci = "../data/pileup/GL-WT-24hr-3.pile_up3_s10_u.shuffled.bed", bin_size = 10, verbose = F) + scale_y_continuous(limits=c(0,45))

fn <- build_filenames(data_path,replicates[c(7,8,9,16,17,18)])
p2 <- plot_bw_profile(bwfiles = fn$files,labels = fn$labels,mode = "center",upstream = 1000, downstream = 1000, loci = "../data/pileup/GL-WT-24hr-3.pile_up3_s10_u.uniq.bed", bin_size = 10, verbose = F) + scale_y_continuous(limits=c(0,45))

embed_plot(p1,"CAD-peaks")

Saving 7 x 7 in image

Download CAD-peaks.pdf

embed_plot(p1,"CAD-peaks.rep")

Download CAD-peaks.rep.pdf

embed_plot(p1,"CAD-ctrl")

Download CAD-ctrl.pdf

ggarrange(p1, p2, p3, ncol = 3, nrow = 1)

Heatmaps of the same regions

fn <- build_filenames(data_path,samples)
p1 <- plot_bw_heatmap(bwfile = fn$files[1], mode = "center",upstream = 1000, downstream = 1000, loci = "../data/pileup/GL-WT-24hr-3.pile_up3_s10_u.uniq.bed", bin_size = 10, verbose = F,zmax = 50)

p2 <- plot_bw_heatmap(bwfile = fn$files[2], mode = "center",upstream = 1000, downstream = 1000, loci = "../data/pileup/GL-WT-24hr-3.pile_up3_s10_u.uniq.bed", bin_size = 10, verbose = F,zmax = 50)

p3 <- plot_bw_heatmap(bwfile = fn$files[3], mode = "center",upstream = 1000, downstream = 1000, loci = "../data/pileup/GL-WT-24hr-3.pile_up3_s10_u.uniq.bed", bin_size = 10, verbose = F,zmax = 50)

p4 <- plot_bw_heatmap(bwfile = fn$files[4], mode = "center",upstream = 1000, downstream = 1000, loci = "../data/pileup/GL-WT-24hr-3.pile_up3_s10_u.uniq.bed", bin_size = 10, verbose = F,zmax = 50)

p5 <- plot_bw_heatmap(bwfile = fn$files[5], mode = "center",upstream = 1000, downstream = 1000, loci = "../data/pileup/GL-WT-24hr-3.pile_up3_s10_u.uniq.bed", bin_size = 10, verbose = F,zmax = 50)

p6 <- plot_bw_heatmap(bwfile = fn$files[6], mode = "center",upstream = 1000, downstream = 1000, loci = "../data/pileup/GL-WT-24hr-3.pile_up3_s10_u.uniq.bed", bin_size = 10, verbose = F,zmax = 50)

ggarrange(p1, p2, p3, p4, p5, p6, ncol = 6, nrow = 1)

CTCF sites - spacing relative to nucleosomes and linker histone

around CTCF sites, the GLOE-seq signal looked very periodic, which suggested that there could be a relationship to nucleosome occupancy. Indeed, the pattern is consistent with CAD cutting left and right of linker histone / in-between linker histone and nucleosome. Note that the same periodicity is also seen in untreated and CAD KO conditions, albeit much less pronounced, suggesting that CAD-independent SSBs are also occuring more often in-between nucleosomes.

DNase hypersensitive sites - spacing relative to nucleosomes and linker histone

With the removal of artifactual regions, we now also can see a neatly spaced pattern at DNAse hypersensitive sites.

p4<- plot_bw_profile(bwfiles = c(plus$files[3],minus$files[3]),labels = c("WT 24h +","WT 24h -"),mode = "center",upstream = 1000, downstream = 1000, loci = "../genome/HCT116_DHS.hg19.bed", bin_size = 10, verbose = F, remove_top = 0.0001)

p4<- plot_bw_profile(bwfiles = c(plus$files[3],minus$files[3]),labels = c("WT 24h +","WT 24h -"),mode = "center",upstream = 1000, downstream = 1000, loci = "../genome/HCT116_DHS.hg19.bed", bin_size = 10, verbose = F, remove_top = 0.0001)

p5 <- plot_bw_profile(bwfiles = plus$files[c(1,3,6)],labels = plus$labels[c(1,3,6)],mode = "center",upstream = 1000, downstream = 1000, loci = "../genome/HCT116_DHS.hg19.bed", bin_size = 10, verbose = F, remove_top = 0.0001)

p6 <- plot_bw_profile(bwfiles = minus$files[c(1,3,6)],labels = minus$labels[c(1,3,6)],mode = "center",upstream = 1000, downstream = 1000, loci = "../genome/HCT116_DHS.hg19.bed", bin_size = 10, verbose = F, remove_top = 0.0001)

ggarrange(p1, p2, p3, p4, p5, p6, ncol = 3, nrow = 2)